R-Spatial Notebook

Spatial data everywhere!

Adding geo dimension

Geographical data is everywhere. In Statoil we use maps to keep track of the licenses, well/installation locations, we use maps to outline the prospects and plan oil/gas pipelines. Our internal ArcGIS database contains a lot of spatial data produced internally (shapefiles for basins, licenses, prospects/fields), as well as procured from external sources (IHS, Tellus, NPD, etc).

However you may still find yourself feeling you are “starting from scratch”. The truth is, you probably have geo data even if you don’t have column for longitude and latitude. These days, google knows everything, and if google knows it, you probably have coordinates.

NB! Google allows only 2500 requests per day through their API free of charge

library(ggmap)
library(dplyr)
statoil_locations <- c("statoil fornebu", "statoil forus", "statoil kontor stjoerdal", "statoil harstad", "statoil tjeldbergodden", "statoil mongstad", "statoil kaarstoe", "statoil melkoeya") 
stl_loc_df <- statoil_locations %>% 
  geocode() %>% bind_cols(location=statoil_locations)

Information from URL : http://maps.googleapis.com/maps/api/geocode/json?address=statoil%20fornebu&sensor=false
Information from URL : http://maps.googleapis.com/maps/api/geocode/json?address=statoil%20forus&sensor=false
.Information from URL : http://maps.googleapis.com/maps/api/geocode/json?address=statoil%20kontor%20stjoerdal&sensor=false
.Information from URL : http://maps.googleapis.com/maps/api/geocode/json?address=statoil%20harstad&sensor=false
geocode failed with status OVER_QUERY_LIMIT, location = "statoil harstad".Information from URL : http://maps.googleapis.com/maps/api/geocode/json?address=statoil%20tjeldbergodden&sensor=false
.Information from URL : http://maps.googleapis.com/maps/api/geocode/json?address=statoil%20mongstad&sensor=false
geocode failed with status OVER_QUERY_LIMIT, location = "statoil mongstad".Information from URL : http://maps.googleapis.com/maps/api/geocode/json?address=statoil%20kaarstoe&sensor=false
.Information from URL : http://maps.googleapis.com/maps/api/geocode/json?address=statoil%20melkoeya&sensor=false

stl_loc_df

Let’s see if we can plot these:

library(leaflet)
p <- leaflet(stl_loc_df) %>% 
  addTiles() %>% 
  addMarkers(popup = ~location) 

Assuming 'lon' and 'lat' are longitude and latitude, respectively
Data contains 2 rows with either missing or invalid lat/lon values and will be ignored

p

Not bad! Lets see if we can get some shapefiles to play with:

library(raster)

Loading required package: sp

Attaching package: 㤼㸱raster㤼㸲

The following object is masked from 㤼㸱package:dplyr㤼㸲:

    select

norway_fylke <- getData('GADM',
  country = 'NO', 
  level = 1) 
p %>% 
  addPolygons(data = norway_fylke)

Of course you can read shapefiles (e.g. Statoil/Tellus basin outlines) and data frames with pre-populated coordinates (e.g. well locations in Woodmac Upstream Data Tool).

Manipulating spatial data

Now that we have spatial data, how do we manipulate it? For example, how do we find out, in which fylke is each of our offices located?

First we will need to transform our data into the “tidy” format (encode geo-data as as simple feature and include it into the “rectangular” dataframe as special “geometry” column). Once the data is in the “tidy” format, we can easily make “spatial join” using familiar “tidyverse” tools.

Please, review awesome vignettes to sf package

library(sf)
stl_loc_geo <- stl_loc_df %>% na.omit() %>% st_as_sf(coords=c("lon", "lat"), crs=4326)
stl_loc_geo

Simple feature collection with 6 features and 1 field
geometry type:  POINT
dimension:      XY
bbox:           xmin: 5.523974 ymin: 58.89266 xmax: 23.60045 ymax: 70.68937
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs
                  location                      geometry
1          statoil fornebu POINT (10.6300284 59.8958549)
2            statoil forus  POINT (5.7190168 58.8926611)
3 statoil kontor stjoerdal POINT (10.9074387 63.4701076)
5   statoil tjeldbergodden  POINT (8.6937046 63.4084316)
7         statoil kaarstoe  POINT (5.5239741 59.2789269)
8         statoil melkoeya  POINT (23.600448 70.6893659)

norway_fylke_geo <- norway_fylke %>% st_as_sf(crs=4326)
stl_loc_geo %>% 
  st_join(norway_fylke_geo, join=st_within) %>% 
  as_tibble() %>% 
  dplyr::select(location, NAME_1)

although coordinates are longitude/latitude, it is assumed that they are planar

What if we want to take geo data from the user, say I want to mark up my favorite hiking areas? (or indicate boundaries of a new prospect or a prospective area to be allocated in the new licensing round).

hiking_areas_geo

Simple feature collection with 3 features and 2 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: 8.9648 ymin: 59.9715 xmax: 10.7364 ymax: 63.4352
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs
  X_leaflet_id feature_type                       geometry
1           81      polygon POLYGON ((10.676 60.0127, 1...
2           94      polygon POLYGON ((9.1406 60.8449, 8...
3          107      polygon POLYGON ((10.2571 63.4302, ...

Ok. Where are those located? Or rather, give me the list of hikeable fylkes (fylke that contain my hiking spots). Note that I marked the areas without looking at the fylke borders so I am interested in a generic join of polygons that may be intersecting.

hiking_areas_geo %>% 
  st_join(norway_fylke_geo) %>% 
  as_tibble() %>% 
  dplyr::select(NAME_1) %>% 
  distinct

although coordinates are longitude/latitude, it is assumed that they are planar

Finally, last summer, my friend Adrian took me fishing outside of Linnesøya in Trøndelag. We had fancy fishing equipment and a sonar device in the boat, so we were able to get a lot of fish, but I always wondered how deep was the water in that area and whether I could go with my regual fishing rod to catch that steinbit that cut off my lure.

cent

         X        Y
1 10.24712 64.41337

so what was the waterdepth in the area we were fishing (using centroid location)? How far was it from shore? What does the seabed look like? First, we fetch bathymetry data from free online service to explore seabed in this location